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■ Abstract 

o 

I We have employed the multiple image method to compute the interparticle 

force for a polydisperse electrorheological (ER) fluid in which the suspended 

^ I particles can have various sizes and different permittivites. The point-dipole 

(PD) approximation being routinely adopted in computer simulation of ER 

^ ^ fluids is shown to err considerably when the particles approach and finally 

\^ ' touch due to multipolar interactions. The PD approximation becomes even 

(N ■ 

I worse when the dielectric contrast between the particles and the host medium 

is large. From the results, we show that the dipole-induced-dipole (DID) 
I model yields very good agreements with the multiple image results for a wide 

, range of dielectric contrasts and polydispersity. As an illustration, we have 

Q I employed the DID model to simulate the athermal aggregation of particles 

in ER fluids both in uniaxial and rotating fields. We find that the aggrega- 
^ ' tion time is significantly reduced. The DID model accounts for multipolar 

interaction partially and is simple to use in computer simulation of ER fluids. 

PACS Number(s): 83.80.Gv, 82.70.Dd, 41.20.-q 
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I. INTRODUCTION 



The prediction of the yield stress for electrorheological (ER) fluids is the main concern in 
theoretical investigations of ER fluids. Early studies failed to derive the experimental yield 
stress data [|l| because these studies were almost based on a point-dipole approximation 
The point-dipole approximation is routinely adopted in computer simulation because 
it is simple and easy to use. Since many-body and multipolar interactions between particles 
have been neglected, the strength of ER effects predicted by this model is of an order lower 
than the experimental results. Hence, substantial effort has been made to sort out more 
accurate models. 

Klingenberg and coworkers developed empirical force expression from numerical solution 
of Laplace's equation 0]. Davis used the finite-element method [^]. Clercx and Bossis 
developed a full multipolar treatment to account for multipolar polarizability of spheres 
up to 1,000 multipolar orders 0]. Yu and coworkers developed an integral equation method 
which avoids the match of complicated boundary conditions on each interface of the particles 
and is applicable to nonspherical particles and multimedia 0. Although the above methods 
are accurate, they are relatively complicated to use in dynamic simulation of ER fluids. 
Alternative models have been developed to circumvent the problem: the coupled-dipole 
model and the dipole-induced-dipole model [0, which take care of mutual polarization 
effects when the particles approach and finally touch. 

The DID model accounts for multipolar interactions partially and is simple to use in 
computer simulation of ER fluids P]. As an illustration, we employed the DID model to 
simulate the athermal aggregation of particles in ER fluids both in uniaxial and rotating 
fields. We find that the aggregation time is significantly reduced. In the next section, we 
review the multiple image method and establish the dipole-induced dipole (DID) model. In 
section III, we apply the DID model to the computer simulation of ER fluids in a uniaxial 
field. In section IV, we extend the simulation to athermal aggregation in rotating fields. 
Discussion on our results will be given. 
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II. IMPROVED MULTIPLE IMAGE METHOD 



Here we briefly review the multiple images method and extend the method slightly 
to handle different dielectric constants. Consider a pair of dielectric spheres, of radii a and 
b, dielectric constants ei and e[ respectively, separated by a distance r. The spheres are 
embedded in a host medium of dielectric constant €2- Upon the application of an electric 
field Eo, the induced dipole moment inside the spheres are respectively given by (SI units): 



(1) 



where the dipolar factors /?, j3' are given by: 
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ei + 2e2 ' " e'l + 2e2 
From the multiple image method 0, the total dipole moment inside sphere a is: 
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where the subscripts T (L) denote a transverse (longitudinal) field, i.e., the applied field is 
perpendicular (parallel) to the line joining the centers of the spheres. Similar expressions 
for the total dipole moment inside sphere b can be obtained by interchanging a and b, as 
well as P and The parameter a satisfies: 
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In Ref. 0, we checked the validity of these expressions by comparing with the integral 
equation method. We showed that these expression are valid at high contrast. Our improved 
expressions will be shown to be good at low contrast as well (see below). 
The force between the spheres is given by ||10|| : 
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For mono disperse ER fluids {a = b, P = P' and Pa = Pb = Po), Klingenberg defined an 
empirical force expression [^]: 

F 



{2K\\ cos' e-K^ sin^ e)r + Kt sin 299, (6) 

FpD 

being normalized to the point-dipole force Fpo = — Spg/r^, where K\\,K± and Kr (all 

tending to unity at large separations) are three force functions being determined from the 

numerical solution of Laplace's equation. The Klingenberg's force functions can be shown 

to relate to our multiple image moments as follow (here a = b, P = P' and pa = pt)'- 

IdpL dpT l,^ . . , . 

^11 = 2^' ^^ = ~^' ^r = -(pr-PL), (7) 

where pi = Pl/FpdEq and pt = Pt/Fp^Eq are the reduced multiple image moments of 
each sphere. We computed the numerical values of these force functions separately by the 
approximant of Table I of the second reference of Ref . and by Eq. (0) . 

In Fig.l, we plot the multiple image results and the Klingenberg's empirical expressions. 
We show results for the perfectly conducting limit {P = 1) only. For convenience, we 
define the reduced separation a = r/{a + b). For reduced separation a > 1.1, simple 
analytic expressions were adopted by Klingenberg. As evident from Fig.l, the agreement 
with the multiple image results is impressive at large reduced separation a > 1.5, for all 
three empirical force functions. However, significant deviations occur for a < 1.5, especially 
for K\\. For a < 1.1, alternative empirical expressions were adopted by Klingenberg. For 
K±, the agreement is impressive, although there are deviations for the other two functions. 
From the comparison, we would say that reasonable agreements have been obtained. Thus, 
we are confident that the multiple image expressions give reliable results. 

The analytic multiple image results can be used to compare among the various models 
according to how many terms are retained in the multiple image expressions: (a) point- 
dipole (PD) model: n = 1 term only, (b) dipole-induced-dipole (DID) model: = 1 to 
n = 2 terms only, and (c) mult ipole- induced- dip ole (MID) model: n = 1 to n = oo terms. 

In a previous work , we examine the case of different size but equal dielectric constant 
{P = P') only. Here we focus on the case a = b and study the effect of different dielectric 



constants. In Fig. 2, we plot the interparticle force in the longitudinal field case against the 
reduced separation a between the spheres for (a) (3 = 9/11 (ei/e2 = 10) and (b) (3 = 1/3 
(ei/e2 = 2) and various (3' / (3 ratios. At low contrast, the DID model almost coincides with 
the MID results. In contrast, the PD model exhibits significant deviations. It is evident 
that the DID model generally gives better results than PD for all polydispersity. 



III. COMPUTER SIMULATION IN THE DIPOLE-INDUCED-DIPOLE MODEL 

The multiple image expressions [Eqs.(3)-(4)] allows us to calculate the correction factor 
defined as the ratio between the DID and PD forces: 

^O) - (^2 _ ^2)4 (^2 _ ^2)4 + (^2 _ ^2 _ ^2)4 ' ^ ) 
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where F^^ = SpaoPfeo/^^, -^pd = -^PaoPwIr"^ and F^^ = -'^PaoPbo/r^ are the point-dipole 
forces for the transverse, longitudinal and F cases respectively. These correction factors 
can be readily calculated in computer simulation of polydisperse ER fluids. The results 
show that the DID force deviates significantly from the PD force at high contrast when (3 
and 13' approach unity. The dipole induced interaction will generally decrease (increase) the 
magnitude of the transverse (longitudinal) interparticle force with respect to the point-dipole 
limit. 

For simplicity, we consider the case of two equal spheres of radius a, initially at rest and 
at a separation da. An electric field is applied along the line joining the centers of the sphere. 
The equation of motion is given by: 

f-^-P* (11) 

where z is the displacement of one sphere from the center of mass. The separation between 
the two spheres is thus d — 2z and the initial condition is d — dQ &X, t — Eq.(ll) 



is a dimensionless equation. We have chosen the following natural scales to define the 
dimensionless variables: 

length ~ a, time = iq ~ — 7; — , lorce = i^o ~ : , 

where Eq is the field strength, m is the masss, rjc is the coefficient of viscosity. Using typical 
parameters, we find to is of the order milliseconds. We have followed Klingenberg to 
ignore the inertial effect, captured by the parameter G. 

The neglect of G can be justified as follows. For values common to ER suspension: r/c ~ 0.1 
Pa s, m ^ 8 X 10^^^ kg, a ~ 5 x lO^^m 0], The inertial term G is of the order 10~^. We 
also neglect the thermal motion of the particles which is a valid assumption at high fields. 
We should remark that the initial separation is related to the volume fraction 0, defined 
as the ratio of the volume of the sphere to that of the cube which contains the sphere , 
i.e. = l^sphere/Kube, and 

2a V60y 

For the PD approximation, Eq.(ll) admits an analytic solution: 
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We integrate the equation of motion by the 4th order Runge-Kutta algorithm, with time 
steps 5t = 0.01 and 0.001 for small and large volume fractions respectively. In Fig. 3(a), we 
plot the displacement d/a versus time graph for the PD case and find excellent agreement 
between analytic and numerical results. 

For the DID model, we have to integrate the equation of motion numerically. In Fig.4, we 
plot the displacement d/a versus time graph for the aggregation of two spheres in uniaxial 
fields. At small volume fractions, i.e., when the initial separation is large, the time for 
aggregation is large and the DID results deviate slightly from the PD results. However, at 
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large volume fractions, the DID results are significantly smaller than the PD calculations. 
The effect becomes even more pronounced at large j3. 

In Fig. 5 (a), we plot the ratio of aggregation time of the DID to PD cases. The results 
showed clearly that the aggregation time has been significantly reduced when mutual polar- 
ization effects are considered. The reduction in aggregation time becomes even pronounced 
for small initial separations. 



IV. ATHERMAL AGGREGATION IN THE ROTATING FIELD 



Recently, Martin and coworkers |jTT| demonstrated athermal aggregation with the rotat- 



ing field. When a rotating field is applied in the x-y plane at a sufficiently high frequency 
that particles do not move much in one period, an average attractive dipolar interaction is 
created. The result of this is the formation of plates in the x-y plane. Consider a rotating 
field applied in the x-y plane: E^j. = Eq cos ut, Ey = Eosinut. The dimensionless equation 
of motion for the two sphere case becomes: 

—- = F\\ cos^ Lut + F± sin^ tut, = — Fr sin 2tt;t, (13) 
dt ' dt ' ^ ' 

where (x, y) is the displacement of one sphere from the center of mass. For large u;, we may 

safely neglect the y component of the motion. In the PD approximation, F\\ = —Qp^/r'^ and 

Fj_ = Sp^/r^, we find the analytic result: 
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The separation between the two spheres is just d = 2x, with the initial separation d = do a.t 
t = 0. In the rotating field case, we also integrate the equation of motion by the 4th order 
Runge-Kutta algorithm, but with 6t = l/{Auj) and l/(40ti;) as the time steps. Note that 
6t = l/{4:Lj) should be the largest time step which can be used because we must at least 
go through a cycle consisting of the transverse and longitudinal field cases. The oscillating 
effect of a rotating field is less observable when the time step is smaller than this maximum 
value. 



In Fig. 3(b), we plot the displacement versus time graph for the PD case in rotating field 
and find an excellent agreement between analytic and numerical results. It is evident that 
the aggregation time is 4 times of that of the uniaxial field case. In fact, Eq.(14) reduces to 
Eq.(12) as a; — > 0. However, at large a;, Eq.(14) becomes: 



That is, in the PD approximation, the time average force becomes 1/4 of that of the uniaxial 
field case. It is because the two dipole moments spend equal times in the transverse and 
longitudinal orientations, while F\\ = —2F± in the PD case, leading to an overall attractive 
force that is 1/4 of the force of the uniaxial field case. When the multiple image force 
is included, we expect that the magnitude of Fy increases while that of F± decreases and 
we expect an even larger attractive force when the spheres approach. In this case, the 
aggregation time must be reduced even more significantly. 

In Fig. 5(b), we plot the ratio of aggregation time of the DID to PD cases for /3 = 1 and 
several uj. The a; = curve is just for the uniaxial field case. The results showed clearly 
that the aggregation time has been significantly reduced when mutual polarization effects 
are considered. The reduction in aggregation time becomes even pronounced for small initial 
separations. It is observed that fluctuations exist when the initial separation between the 
spheres is 2.4a or less. It is because the motion is sensitive to the initial orientation of the 
dipoles when the spheres are too close. 

Similarly, we consider the aggregations of 3 and 4 equal spheres, arranged in a chain, an 
equilateral triangle and a square. For a chain of 3 spheres in a rotating field, the central 
sphere does not move, while the two spheres at both ends move towards the central sphere. 
For 3 spheres in an equilateral triangle, the center of mass (CM) will not move while each 
sphere moves towards the CM, subject to the force of the other two spheres. The same 
situation occurs for 4 spheres in a square, in which each sphere moves towards the CM, 
subject to the force of the other 3 spheres. 

In the PD approximation, we report the analytic results as follows. For 3 spheres in a 
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For 3 spheres in an equilateral triangle, 
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For 4 spheres in a square. 
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In each of the above cases, x is the distance of one sphere from the center of mass. In the 
case of 3 spheres in a chain, the separation between the spheres is same as o? = a;. In the 
case of 3 spheres in an equilateral triangle, the separation between spheres is d — y/Sx. In 
the case of 4 spheres in a square, the separation between spheres is d — \/2x. Again, we 
integrate the equation of motion by the 4th order Runge-Kutta algorithm. We find excellent 
agreements between the analytic and numerical results (not shown here). 

It has been found that the displacement in the y-direction is about 0.5% for a; = 5 and 
a larger lv has been used in the simulation. On the other hand, it is time consuming for 
simulations with a; > 20. It is evident from the displacement-time graph that the results are 
correct. In Fig. 6, the oscillation amplitude is reduced when the rotating frequency increases 
in the simulation. This is consistent with the assumption made in our analytic expressions. 
In Fig. 7, it is observed that fluctuations exist when the initial separation between the spheres 
is 2.4a or less in all three graphs. Again, it is because the motion is sensitive to the orientation 
of the dipoles when the spheres are close. From the simulation, the reduction effects become 
even pronounced for the rotating electric field case than the uniaxial field case. 



DISCUSSION AND CONCLUSION 



Here a few comments on our results are in order. In this work, we studied the aggregation 
time for several particles. We should also examine the morphology of aggregation, due to 
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multiple image forces. In this connection, we can also examine the structural transformation 
by applying the uniaxial and rotating fields simultaneously ||12|| . 

We have done simulation in the monodisperse case. Real ER fluids must be polydisperse 
in nature: the suspending particles can have various sizes or different permittivities. Poly- 
disperse electrorheological (ER) fluids have attracted considerable interest recently because 
the size distribution and dielectric properties of the suspending particles can have significant 



impact on the ER response ||T3|. We should extend the simulation to polydisperse case by 
using the DID model. 

In summary, we have used the multiple image to compute the interparticle force for 
a polydisperse electrorheological fluid. We apply the formalism to a pair of spheres of 
different dielectric constants and calculate the force as a function of the separation. The 
results show that the point-dipole approximation is oversimplified. It errs considerably 
because many-body and multipolar interactions are ignored. The dipole-induced-dipole 
model accounts for multipolar interactions partially and yields overall satisfactory results in 
computer simulation of ER fluids while it is easy to use. 
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FIGURES 

FIG. 1. Comparison of multiple image results with Klingenberg's force functions. 

FIG. 2. Force for longitudinal fielda: (a) P=9/ll, (b) (3=1/3. 

FIG. 3. Comparison between analytic and simulation results in athermal aggregation: (a) 
uniaxial field, (b) rotating field. 

FIG. 4. Diplacement-time graph for athermal aggregation of two spheres in uniaxial field. 

FIG. 5. Reduction of aggregation time for: (a) uniaxial field, (b) rotating field. 

FIG. 6. Diplacement-time graph for athermal aggregation of two spheres in rotating field. 

FIG. 7. Reduction of aggregation time for three and four spheres in rotating field. 
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